vegan包进行微生物群落主坐标分析（PCoA）及ggplot2作图
在R中，可用于进行PCoA分析的R包有很多可供选择，如“vegan”、“ade4”等，这些均是在生态统计中常用的R包。此处作为示例介绍其中的“vegan”包。
首先介绍示例数据。我们此处共有96个16S测序样本，均来自土壤。这96个样本共涉及了4个采样地点（地点A、B、C、D）；2种处理梯度，即添加某化学物质的低浓度处理（low）以及高浓度处理（high）；4个采样时期（时期1、2、3、4）；对于每个采样地的每种处理梯度的每个时期下，各自进行了3个重复，共计4×2×4×3=96。
此处我们希望通过PCoA分析，查看样本间细菌群落组成是否具有显著不同。
作图示例文件、R脚本等，已上传至百度盘，无提取码（若失效请在下方留言）
https://pan.baidu.com/s/1mL92U9DxMZKjqTNRhDW9ZA
示例文件简要
文件“otu_table.txt”为OTU丰度表格，其内容展示如下。
每一列为一个样本，每一行为一种OTU，交叉区域为每种OTU在各样本中的丰度。
每一列为一个样本，每一行为一个样本，交叉区域为样本间的Bray-curtis距离（取值范围0-1，越接近于1表明样本间细菌群落组成差异越大）。
第一列（names）为各样本名称；第二列（site）为各样本的采样地点，即4个采样地点（地点A、B、C、D）；第三列（deal）为2种处理梯度，即添加某化学物质的低浓度处理（low）以及高浓度处理（high）；第四列（time）各样本的4个采样时期（时期1、2、3、4）；第五列（repet）为每个采样地的每种处理梯度的每个时期下各自进行的3个重复（1、2、3）。
使用vegan包进行PCoA排序分析
首先导入数据。我们可选导入原始的OTU丰度表格文件，也可使用已经计算好的样本距离矩阵文件，同时导入样本分组文件。
#OTU 丰度表otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)otu <- data.frame(t(otu))#或者现有的距离矩阵dis <- read.delim('bray.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE) #样本分组文件group <- read.delim('group.txt', sep = '\t', stringsAsFactors = FALSE)然后加载vegan包，并进行PCoA分析。
对于以下命令的更多参数，可使用?vegdist()和?cmdscale查看。
library(vegan) #排序（基于 OTU 丰度表）distance <- vegdist(otu, method = 'bray')pcoa <- cmdscale(distance, k = (nrow(otu) - 1), eig = TRUE)#或者（基于现有的距离矩阵）pcoa <- cmdscale(as.dist(dis), k = (nrow(dis) - 1), eig = TRUE)#对于上述计算得到的样本间距离 distance，可转换为矩阵格式后输出，例如输出为 csv 格式
write.csv(as.matrix(distance), ‘distance.csv’, quote = F)
若我们之前导入的是OTU丰度表，则我们需要首先根据OTU的丰度组成，计算样本间距离，然后使用计算好的样本间距离对样本进行PCoA排序。vegdist()用于计算样本间距离，此处使用Bray-curtis距离；cmdscale()用于PCoA排序。eig选择“TRUE”或“FALSE”对结果的影响较大，参数选择需要根据数据分布、所关注的生物学意义等因素慎重考虑。
若我们之前导入的是现有的距离矩阵，则我们可直接基于先有距离对样本进行PCoA排序。首先使用as.dist ()转化读入的矩阵，并使用cmdscale()进行PCoA排序。
此处展示了两种过程，若两种过程中所使用的“距离类型”是一致的，则最后所得结果也是相同的。如此处展示的两种方式均使用了Bray-curtis距离，因此排序结果完全一致。
可以简要地查看结果。
#使用vegan内置命令ordiplot()简要做图展示ordiplot(scores(pcoa)[ ,c(1, 2)], type = ‘t’)#或者查看排序简要summary(pcoa)vegan内置命令ordiplot()方便我们直接查看排序结果。若结果可观，我们可以考虑将排序结果中的各项指标提取出，绘制效果更好的排序图（例如借助ggplot2）。
summary(pcoa)的打印结果中，给出了排序结果中的各项重要指标。
主要关注两个重要指标，eig记录了PCoA排序结果中，主要排序轴的特征值（再除以特征值总和就是各轴的解释量）；points记录了各样本在各排序轴中的坐标值。
#查看主要排序轴的特征值和各样本在各排序轴中的坐标值pcoaKaTeX parse error: Expected 'EOF', got '&' at position 10: eigpoint&̲nbsp;<- …point)#可将样本坐标转化为数据框后导出，例如导出为 csv 格式write.csv(point, ‘pcoa.sample.csv’)
此时我们得到了各样本的PCoA排序图。
我们还可使用vegan包中的命令wascores()，得到各OTU的排序坐标。wascores()可以以多度加权平均方式将各OTU被动投影到样本的PCoA排序图内，可使用?wascores()查看详情。因OTU数据量较大，因此在这里只展示前两个排序轴。
#可使用 wascores() 计算物种坐标species <- wascores(pcoaKaTeX parse error: Expected 'EOF', got '&' at position 14: points[,1:2],&̲nbsp;otu) …eig)[1:2] / sum(pcoaKaTeX parse error: Expected 'EOF', got '&' at position 6: eig)&̲nbsp;#提取样本点坐标（…point})[1:2]sample_siteKaTeX parse error: Expected 'EOF', got '&' at position 6: names&̲nbsp;<- …site <- factor(sample_siteKaTeX parse error: Expected 'EOF', got '&' at position 6: site,&̲nbsp;levels&nbs…deal <- factor(sample_siteKaTeX parse error: Expected 'EOF', got '&' at position 6: deal,&̲nbsp;levels&nbs…time <- factor(sample_site$time, levels = c(‘1’, ‘2’, ‘3’, ‘4’))library(plyr)group_border <- ddply(sample_site, ‘site’, function(df) df[chull(df[[2]], df[[3]]), ])#注：group_border 作为下文 geom_polygon() 的做图数据使用然后使用ggplot2进行PCoA排序图绘制。
此处分组较多，因此在本示例中，考虑使用多边形区域展示不同采样来源（A、B、C、D四个地点，绘制方法可参考前述博文），使用两种形状区分2种梯度的处理（low、high），使用渐变颜色区分4个采样时期（1、2、3、4四个时期）。
library(ggplot2)pcoa_plot <- ggplot(sample_site, aes(PCoA1, PCoA2, group = site)) +theme(panel.grid = element_line(color = ‘gray’, linetype = 2, size = 0.1), panel.background = element_rect(color = ‘black’, fill = ‘transparent’), legend.key = element_rect(fill = ‘transparent’)) + #去掉背景框geom_vline(xintercept = 0, color = ‘gray’, size = 0.4) + geom_hline(yintercept = 0, color = ‘gray’, size = 0.4) +geom_polygon(data = group_border, aes(fill = site)) + #绘制多边形区域geom_point(aes(color = time, shape = deal), size = 1.5, alpha = 0.8) + #可在这里修改点的透明度、大小scale_shape_manual(values = c(17, 16)) + #可在这里修改点的形状scale_color_manual(values = c(‘yellow’, ‘orange’, ‘red’, ‘red4’)) + #可在这里修改点的颜色scale_fill_manual(values = c(’#C673FF2E’, ‘#73D5FF2E’, ‘#49C35A2E’, ‘#FF985C2E’)) + #可在这里修改区块的颜色guides(fill = guide_legend(order = 1), shape = guide_legend(order = 2), color = guide_legend(order = 3)) + #设置图例展示顺序labs(x = paste('PCoA axis1: ', round(100 pcoa_eig[1], 2), ‘%’), y = paste('PCoA axis2: ', round(100 pcoa_eig[2], 2), ‘%’)) +
#可通过修改下面四句中的点坐标、大小、颜色等，修改“A、B、C、D”标签annotate(‘text’, label = ‘A’, x = -0.31, y = -0.15, size = 5, colour = ‘#C673FF’) +annotate(‘text’, label = ‘B’, x = -0.1, y = 0.3, size = 5, colour = ‘#73D5FF’) +annotate(‘text’, label = ‘C’, x = 0.1, y = 0.15, size = 5, colour = ‘#49C35A’) +annotate(‘text’, label = ‘D’, x = 0.35, y = 0, size = 5, colour = ‘#FF985C’)ggsave(‘PCoA.png’, pcoa_plot, width = 6, height = 5)ggplot2的各细节不再详细说明，最终输出结果如下。
后期可使用AI、PS等工具进行加工。
参考文献
DanielBorcard, FranoisGillet, PierreLegendre, et al. 数量生态学:R语言的应用（赖江山 译）. 高等教育出版社, 2014.